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ABSTRACT 

We use the Fokker-Planck equation and model the dispersive dynamics of soUd particles in annu- 
lar protoplanctary disks whose gas component is more massive than the particle phase. We model 
particle-gas interactions as hard sphere collisions, determine the functional form of diffusion coeffi- 
cients, and show the existence of global unstable modes in the particle phase. These modes have spiral 
patterns with the azimuthal wavenumber m ~ I and rotate slowly. We show that in ring-shaped disks, 
solid particles subject to gas drag never fall into the central star and instabilities occur for particles 
of all sizes. Therefore, planetesimals and planetary cores can be efficiently produced through the 
accumulation of smaller objects near the peaks of unstable density waves. 

Subject headings: methods: numerical, hydrodynamics, instabilities, planets and satellites: formation, 
planetary systems: protoplanctary disks 
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1. INTRODUCTION 

Protoplanctary disks are multi-phase environments 
composed of solid particles and molecular gas. The mo- 
tion of particles is mainly governed by the gravitational 
forces of the disk material and the central star. Gas 
molecules feel the pressure gradient as well: for a poly- 
tropic isothermal gas, whose density profile monotoni- 
cally increases towards the central star, the pressure gra- 
dient Vp is always negative. This yields a sub-Keplerian 
circular velocity and generates headwind on solid par- 
ticles that move on Keplerian orbits. Solid particles 
are thus expected to inspiral towards the central star. 
Although adhesive and electrostatic for ces can enhance 
the c lustering of micron-sized particles (jBlum fc WiirmI 
centimeter- to meter-size particles seem to inspiral 
towards (and fall into) the central star sooner than the 
time scale needed for assembling planetary cores. There- 
fore, several col lective processes like streaming instability 
ijYoudin &: Goodma n .2005), turb ulent vortices induced 
by Kelvin-Helmhol tz instability (jJohansen et al.l 120061 : 
iBai fc Stondl2010airbh and m agnetorotational instability 
pohansen et al.l 120071 120 lit) have been proposed to be 
responsible for the formation of km-scale planetesimals. 

The inspiraling motion of particles, however, does not 
globally occur in disks with non-monotonic density pro- 
files, which are likely to form through a combination of 
viscous ac cretion and photoevapq ration by the central 
star (e.g.. iMatsuvama et al.l 120031 ). A ring-like disk is 
the simplest system with non-monotonic density profile. 
An interesting property of such systems is that Vp be- 
comes positive in regions where the density profile is ris- 
ing, and gas molecules move with super-Keplerian veloc- 
ities. Consequently, solid particles that approximately 
move on Keplerian orbits are accelerated by gas and mi- 
grate outwards. This means that solid particles do not 
necessarily fall into the central star and may instead mi- 
grate into regions where both the head and tail winds 
are minimized. Possible instabilities triggered by such 
migrations can have remarkable implications for the for- 
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mation of planetesimals and planetary cores. 

In this paper , we generalize the analysis of 
iJalali fc Tremaind ()2012l hereafter JT12) to disks that 
include a locally isothermal gas component, and search 
for global instabilities in the particle phase. Our disks are 
self-gravitating and the particle phase has non-zero radial 
velocity dispersion. The dynamics of particles is modeled 
by the Fokker-Planck equation, and the gas component 
is assumed to be in a steady-state rotation around the 
central star. For simplicity, we confine our study to disks 
with Eg/Ep ^ 1 where Sg and Ep are the surface densi- 
ties of the gas and particle phases, respectively. By this 
assumption, the gas component does not respond to the 
disturbances of the particle phase. We neglect collisions 
between solid particles, but those between gas molecules 
and solid particles are taken into account. 

2. DISK MODEL 

We assume a two-phase medium consisting of solid par- 
ticles and gas molecules, and refer to them by subscripts 
p and g, respectively. Solid particles are assumed to be 
monodisperse hard spheres of mass nip and radius Tp. 
The gas phase has the molecular mass mg, and the av- 
erage radius of gas molecules is Tg. We work with ini- 
tially axisymmetric disks whose gas component has not 
streaming motion in the radial direction, but particles 
can move on eccentric orbits. For both the particle and 
gas phases we use the annular ring model of JT12 whose 
dimcnsionless surface density is 



So(i?) = 



3yU 



i?2 



in (1 + i?2)5/2 



R 



r 
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where Afd = A/p + Mg is the total disk mass. Mi, is the 
mass of central star, Mp is the mass of particles, Mg is the 
mass of gas component, 6 is a length scale, and r is the 
radial distance to the central star. Figure [T] shows the ra- 
dial profile of Eq/m that peaks at i? = VG/B. We suppose 
that the surface densities of the particle and gas phases 
are proportional to Eg so that Ey(i?) = AfyEo(i?) with 
= p,g and — My/M^. Defining G as the gravita- 
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Figure 1. Density profile of ring-shaped disks. 



tion constant, the dimensionless total gravitational po- 
tential field arising from Sq and the central star becomes 



VoiR) = 



1 

R 



H 1 + 2R^ 

2 (l+i?2)3/2' 



(2) 



where $0 is the actual potential. 

We work with initially axisymmetric disks whose gas 
component has not streaming motion in the radial direc- 
tion. If the gas is locally isothermal, its dimensionless 
pressure will be determined from 



P 



/ b 



knT„ 



(3) 



where is the normalized sound speed, fee is the Boltz- 
mann constant, and Tg is the absolute temperature of 
the gas. In passive disks heated by the stellar radiation 
(not by accretion), and sufficiently far from the central 
star, the radial profile of Tg is given by ()Armitagell201Cll 
§2.4) 
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where is the effective temperature of the central star 
and Vi, is its radius. The mean thermal speed of gas 
molecules is related to the sound speed as u^j^ — {8/tt)c^. 
We therefore find 



Uth 



1/2 



\27rb) ' 
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In the solar system, with b set to Jupiter's orbital semi- 
major axis, and nig being the molecular mass of H2, one 
obtains a « 0.022. For Fomalhaut's ring with b « 140AU 
and for the inner ring of AB Aurigae with b « 40 AU 
(|Hashimoto et al.l[201ll ). we find a « 0.036 and 0.026, 
respectively. 

In the absence of collisions between solid particles and 
gas molecules, the velocity of solid particles on circular 
orbits is determined from 
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and the radial momentum equation for the gas becomes 



R 



dR ^ dR' 



from which we obtain the circular velocity of gas: 



Stt (1 -3i?2)^ 2 
24 (1 + 



(7) 
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The condition ^> implies ji < 5^^^, which is satisfied 
by protoplanetary disks. It is seen that for R^ < 1/3, the 
gas streaming velocity exceeds the speed of particles and 
generates tailwind on them. This is a remarkable and 
largely overlooked feature of ring-shaped disks, and has 
interesting implications for the dynamics of solid parti- 
cles in protoplanetary disks: while particles experience a 
resistive headwind for i?^ > 1/3 and inspiral towards the 
central star, they are accelerated when R^ < 1/3 and mi- 
grate outwards. Inward and outward migrating particles 
will then be accumulated around R^ « 1/3, but this is 
not the whole story. In next sections we show that global 
exponentially growing instabilities accompany such ra- 
dial migrations. 

3. EVOLUTIONARY DYNAMICS OF SOLID PARTICLES 

The dynamics of particles is described by the 
phase space distribution function (DF) f{x,v,t) = 
mpj\fp{x,v,t) where J\fp is the number density of par- 
ticles. The vectors x = {xi,X2) and v — (ui,U2) (both 
in Cartesian coordinates) are the position and velocity 
vectors of particles in the disk plane, and t is the time. 
We utilize the DF (jJalah fc Tremainell2012i §3) 



(9) 



to model the initial distribution of particles, before turn- 
ing on the collisions between solid particles and gas 
molecules. Here L = \x x v\ and E = ■ v + Vq are, 
respectively, the orbital angular momentum and energy 
of particles per unit mass. K is a positive integer that 
controls the mean eccentricity e of the disk. The function 
gK{E) takes physical (positive) values ioi K > 2, and e 
decreases as K is increased. In the limit oi K 00 the 
disk becomes cold with all particles moving on circular 
orbits. 

In this study we ignore particle-particle collisions, 
and assume that colliding solid particles and gas 
molecules are hard spheres. The evolution of / 
is therefore expressed by the Fokker-Planck equation 
()Binnev fc Tremaind [20081 §7.4) 



dt ' dxi 
1 92 



+ a, 



dl 

dvi 



d_ 

dvi 



{D[^v,] I) 



2 dvjdv, 



iD[Av,Av,] f) , 



(10) 



where D[/S.Vi] and D[A.Vi/S.Vj] are diffusion coefficients 
and fli arc the components of the acceleration vector. In 
equation (jlOp and throughout the paper a repeated in- 
teger index stands for summation over that index from 
1 to 2. The acceleration vector is computed from a = 
-VT^ = -V[Vb + Vi(i?, (/),*)] where the perturbed po- 
tential Vi is self-consistently calculated from the density 
disturbance Yii[R,<f),t) of particle distribution. Since we 
have assumed Ep ^ Sg, the contribution of the gas com- 
ponent to Vi is neglected. 

For local collisions, the d iffusion coefficient s are e valu- 
ated using the procedure of lRosenbluth et ahl (I1957D . but 
for three dimensional collisions of hard spheres with the 
cross section a = /S^ /4 where P = + r^. Let us define 
7 = nig/ (nip + nig) and denote the absolute velocity of 
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gas molecules by t/. We obtain 

, dh{x,v,t) , , d'^q(x,v,t) , , 



dvi 



dvidvj 



where the potential functions h(x, v,t) and 17(0;, v,t) are 
given by the following integrals 

h^^^^^ f Ms{x,v',t)\v~v'\'^dv', (12) 



9- 



3 
15 



Mg{x,v',t)\v-v'fdv', (13) 



and Mg{x, 1/ ,t) is the number density of gas molecules in 
the phase space. Defining 



(14) 



the elements of the stress tensor are determined as r^j = 
(viVj — ViVj). The macroscopic quantities Sp, Vi and 
Ty are functions of x and t, and the second term on the 
right-hand side of equation (jlOp integrates to zero up 
to the first-order moment equations. CoUisional terms 
involving g{x^v,t) correspond to random motions, and 
contribute to the evolutionary equations of (second- 
order moments of the Fokker-Planck equation). We ne- 
glect them in the present study because the mass ratio 
of gas molecules to solid particles is small, mg/irip <C 1, 
which implies g/h ^ ©(mg/mp). The Fokker-Planck 
equation can therefore be reduced to 



a/ 
at 



d 



d 



dxi 



with 



dV 

dxi 



(15) 



(16) 



Differentiating ([T^ with respect to Vi gives 



TT-lP^ Ng{x,v' ,t)\v ~ v'\{vi ~ v'i) dv' .{n) 



D[Avi 



Evaluating this integral requires the explicit form of 
A/'g, which can be modeled (for example) by a Maxwell- 
Boltzmann distribution. Nonetheless, such details are 
not necessary if we make some further simplifying as- 
sumptions: Let (efl, e^) be unit base vectors in the polar 
coordinates [R^ cj)), and e ^ 1 denote the orbital eccen- 
tricity of particles. The velocities of particles and gas 
molecules will be approximated as 

Vp^ce^ + 0{em), v'-ug^^e^-f 0(z;th), (18) 

where = R~^/^ + 0{^) is the orbital frequency of 
particles. We think of disks with uth < ^p,c, Wg,c over 
0.032 ^ i? ^ 10 (cf. Figure 1) to guarantee the exis- 
tence of bound orbits, and avoid gas dispersal through 
photoevaporation. The mean eccentricity corresponding 
to (O is almost constant over the entire disk space (see 
JT12), anditisgivenby e= ['!r/{4:K + 2)]^/'^ + 0{n). For 
a ^ O(10~^), which corresponds to protoplanetary disks 
around solar-type stars, and for K < 29 experimented in 
JT12, one has 



0{ii), 0.032 <i?< 10. (19) 



From ([5]), (HH) and (HI]) we conclude that the bulk of 
particles satisfy \v — v'\ ^ eRfl, and the integral in (|17p 
is approximated by 



D[Av,]^-^oeRnpg {v,-U,). 



(20) 
(21) 



where pg = Sg(i?)(5(z) is the spatial gas density nor- 
malized to M^/b^, and Ui{x,t) = v[ are the streaming 
velocity components for the gas. 5{z) is the Dirac delta 
function. It is remarked that all velocities are normal- 
ized to [GM^/fo]^/^. To trace the evolution of particles 
in the disk plane, we soften 5{z) using its Poisson kernel 
as 5^{z) = 7r~-^e/(e^ -I- z^) and evaluate D[Avi] at z = 0. 
The small parameter e determines the vertical size of the 
disk. With e = rp/6 the scattering of gas molecules takes 
place in the disk plane and the cross section a becomes 
a line of the length rp -f rg . For near-circular orbits with 
e — >■ 0, one obtains \v— v/\ k, vt\i and D[Avi\ transforms 
to the well-known form of Epstein drag. 

Solving (|15|) in a four-dimensional phase space, with 
particle motions confined to the disk plane, is facilitated 
by utilizing the angle variables w = {wi^W2) and their 
conjugate actions J = (Ji, J2). We follow JT12 and set 
Ji and J2 to the radial action and the orbital angular 
momentum respectively. In the {w, J)-space, the 

motion equations (|16p become 



dU 



Oji OWi 



1,2, (22) 



where T-L = ■ v + V{x, t) is the Hamiltonian function, 
and the functions Fj^ and fu,. are determined using the 
virtual work of nonconservative forces: 



FwiSwi + Fj.SJi = D[Avi] ■ Sxi 



(23) 



with S being the variational operator. After expanding 
D[Avi] and Sxi in Fourier series of angle variables, we ob- 
tain F„. (J, wi) and Fj.{J,'Wi), which are real harmonic 
functions of wi and are smooth in the J-space. We now 
write equation (|15p in the angle-action space: 



(24) 



where [. , .] denotes the Poisson bracket over the (w, J) 
space and the operator I?^ is defined by 

dF^^ dFj^ 



dw. 



4. UNSTABLE MODES 



We seek solutions of the form / = fo{J) + fi{w, J,t) 
for equation ([24)) so that |/i| ^ |/o|. The Hamiltonian 
corresponding to / will become H = Hq + V\ where the 
perturbed potential V\ (self-consistently arising from /i) 
and Ho = \v ■ V + Vb(i?) are expressed in the angle- 
action space (iJalalil [20101 §2). /i is obtained by solving 
the linearized Fokker-Planck equation: 



di 



[/i,Ho] + [/o, V^i] +2?f/i = -2?f/o, (26) 



which is a non-homogenous partial differential equation. 
The particular solution of p6)) is a steady drift of the 
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Table 1 

Pattern speeds and growth rates of slow modes Al 
and A2 versus ^. 





Mode Al 


Mode A2 


UJr X 10-^ 


s X 10^ 


LOr X 10^ 


s X 10^ 





5.172 





4.944 





0.04 


5.172 


1.203 


4.944 


1.631 


0.12 


5.175 


3.613 


4.942 


4.887 


0.20 


5.181 


6.036 


4.939 


8.127 


0.28 


5.190 


8.471 


4.936 


11.34 


0.36 


5.201 


10.91 


4.932 


14.51 


0.50 


5.225 


15.18 


4.927 


19.75 



form fd{J,wi) that satisfies 

[fd,Ho] + [fo,Vd]+VFfd 



--Dpfo, 



(27) 



where the potential Vd corresponds to fd- Wc arc in- 
terested in non-axisymmctric solutions that depend on 
(^1,^2) and t. For a slowly rotating density wave 
Yji{R,(j),t) = Re J fidv with the azimuthal wavenum- 
ber m, we consider unsteady DFs of the form (JT12, 
§6) /i = /i(J) exp[— iwt + im(K;2 — wi)] and compute 
the eigenfrequency a; = oj^ + is and /i through solv- 
ing the homogeneous part of equation p6|) . Here Wr and 
s are, respectively, the pattern speed and growth/decay 
rate of density perturbations. We formulate the eigen- 
value problem following the finite element method of 
[jal ali (201fl) and JT12, and compute the eigenfrequency 
spectrum of m = 1 waves for a model with = 29, 
a = 0.025, ^JL = 0.08, and Mp/Mg = 1/8, and for several 
choices of ^ = £^o/{TTe). The mean eccentricity of parti- 
cle orbits of this model is e = 0.156. Our finite element 
model has = 180 ring elements whose radial nodes are 
located at i?„ = io(4»-2W-4)/(Ar+i) for n = 1, 2, . . . , iV. 
Using this mesh, eigenfrequencies are calculated with a 
relative accuracy < 0.1%. 

In the absence of gas drag, ^ = 0, the spectrum con- 
tains only two stable slow modes, which are labelled 
by Al and A2. The pattern speeds of these modes 



satisfy the inequality oJr > f^p 



0.00469 where 



no 



pr — ^"2 

tide orbits 



rii is the precession frequency of par- 
The orbital frequencies f^i = dHo/dJi 
{i ~ 1,2) are calculated in the unperturbed state. All 
other modes are singular and form a continuum over the 



range ^ipr.min < 



UJr < i^pr,max- By Setting ^ > and 



turning on the gas drag, all singular and long-wavelength 
modes become unstable. The pattern speeds and growth 
rates of modes Al and A2 have been reported in Table 
[1] for several values of ^. Assuming that the mass of 
each particle is computed from nip = (4/3)7rpsrp, with 
Ps being the typical density of rocky material, one finds 
^ (X 1/rp. Therefore, instability- induced accumulation of 
smaller particles happens faster. 

As Table [1] shows, mode A2 grows faster than mode 
Al but it rotates slower. Gas drag has not a notable 
contribution to the pattern speeds of modes and it only 
controls the growth rate so that s linearly increases in 
terms of ^. Figure [2] displays the perturbed density pat- 
terns of stable and unstable modes. It is evident that 
these modes are confined to the region with rising density 
profile. Wc find that the spirality of mode Al increases 
proportional to Comparing the modal content of our 



Mode A1 
^=0 



Mode Al 
^=0.28 



Mode A2 
4 = 



c 



Mode A2 
4 = 0.28 



r 



Figure 2. Mode shapes of global stable (top) and unstable (bot- 
tom) density waves of the particle phase in a protoplanetary disk 
with fi = 0.08, Mp/Mg = 1/8 and e = 0.156. The parameter S. 
is equal to and 0.28 for stable and unstable modes, respectively. 
Contour plots show the positive parts of the perturbed density Si 
at t = 0. Unstable mode Al is a trailing spiral while the wave 
packets of unstable mode A2 form a leading pattern. Filled circle 
shows the location of the central star. 

^ = disks with the results of JT12 (see their Figure 4) 
shows that short-wavelength slow modes have been dis- 
appeared by setting Mp < and the pattern speeds of 
the modes with the longest wavelengths have approached 
to 17^^ jjjg^^. This is very similar to the behavior of disk 
galaxies: increasing the mass of dark matter halo sta- 
bilizes tightly wound spiral modes and only modes with 
the longest w avelengths, especially the bar mode, survive 
(|Jalali|[2007D . The number of slow modes also depends 
on the mean eccentricity of particle orbits. Our numer- 
ical experiments show that mode A2 disappears as its 
pattern speed drops well below ripr,max by increasing e. 

Using equation ([77]) one can show that the particle 
distribution has a secular increase towards R'^ w 1/3 
so that fd{a,e,wi) ~ ©(jywi/rii) ~ 0{rit) where r] = 
£,eail{a)T,g{a), and a is the instantaneous semi-major 
axis of a test particle's orbit. At i? = v^/B where the 
particle density is maximum we have sa 5 x 10~* ^, 
which has the same order of magnitude of s. We con- 
clude that since instabilities grow exponentially as /i ~ 
0(e**), they will overwhelm the particle accumulation 
induced by radial drift. 

5. DISCUSSIONS 

The infall time scale of 0(10^) yr for meter-sized solid 
bodies puts a strong constraint on the fo rmation of plan- 
etesim als from dust grains and pebbles (jWeidenschillingl 
Il977| ) . Recent simulations have shown that the backreac- 
tion of particles on gas can trigger out-of-plane Kelvin- 
Helmholtz instability, which boosts local particle den- 
sity and helps self-gravity to assemble km-size bodies. 
Turbulence, on the other hand, imposes stochastic forces 
on planetesimals, increases their collision frequency and 
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disru pts those with < 10 km radius (jNelson fc Gressel 
120101) . Although the random forcing of planetesimals 
can be suppressed in the presence of a dead zone 
ijGressel et al.l[20TTh . alternative and simpler processes 
may also be involved in the formation of super km-scale 
planetesimals. 

Slow density waves exist in all near-Keplerian, self- 
gravitating rings and can be excited by encounters 
(JT12). Addition of a gas component, however, destabi- 
lizes the particle phase without any external disturbance. 
We showed that the drag-induced infall of solid bodies 
into the central star is not a universal phenomenon and 
the direction of particle migration highly depends on the 
disk structure. There will be no infall if at some re- 
gion the disk density profile, including its solid particle 
and gas components, rises outwards. Therefore, with a 
preserved source of solid particles in the disk, global in- 
stabilities explored in this study will have time to boost 
the particle density to arbitrarily large levels and enhance 
the formation of bigger objects through gravitational col- 
lapse. 

Our results have been obtained for self-gravitating 
disks whose particle phase is constituted from monodis- 
perse hard spheres. Including the size distribution of par- 
ticles and the effect of particle-particle collisions, espe- 
cially catastrophic disruptions, is an interesting problem 
for future studies. Moreover, by assuming Sp <C Eg we 
neglected the backreaction of particles on gas flow. Gen- 
eralization of our study to disks with Sp > Eg requires si- 
multaneous solution of hydrodynamic and Fokker-Planck 



equations. In such circumstances, we anticipate effective 
angular momentum exchange between the particle and 
gas phases, and the gas phase may also become unsta- 
ble. 

I thank Scott Tremaine for his stimulating discussions 
and useful comments. 
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